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Abstract 

In this supplementary material, we recollect, for reader's convenience, the general 
scheme of suggested multiscale model (Sec. [T]), and basic informations about ap- 
proaches used for pilot study: a detailed description of the interaction model (Sec. 
[2]) and dynamical methods used for the dark dynamics step (Sec. [s]) reported pre- 
viously in two preceding studies [D H] . In addition, a detailed description of the 
treatment of radiative processes is also given (Sec. |4]). 
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1 General model scheme 




2 Interaction model 

The intra-cluster interactions are described within an extended diatomics-in- 
molecules (DIM) model with the spin-orbit interaction included. The DIM 
approach was developed by Ellison [3] and later on applied to singly ionized 
rare-gas cluster cations by Kuntz and Valldorf [1]. How the spin-orbit (SO) 
interaction can be included in the DIM model was proposed by Amarouche 
et al. |5]. Note that, in addition to the SO term, further extensions to 
the original DIM approach can also be considered, e.g., the inclusion of the 
leading three-body polarization forces [5]. However, we do not use these 
extensions in the present study as they do not contribute much for small 
clusters [T], and, consequently, we omit them from the following explanation. 

The original DIM approach consists in a) re- writing the electronic hamil- 
tonian as a sum of diatomic and atomic contributions, 

n— 1 n n 

H = E E iijk-in-2)J2Rk, (1) 
j=l k=j+l k=l 
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where n denotes the number of atoms, and b) designing an appropriate 
basis set of wave functions for which the elements of the corresponding 
hamiltonian matrix can be calculated by means of the electronic energies of 
atomic and diatomic fragments. If the SO coupling is not considered, the 
atomic contributions of Eq. [T] are constant and their sum can be identified 
with the zero energy levelj^ The diatomic energies are to be supplied from 
independent sources (usually from ab initio calculations). 

The basis set proposed for an n-atom rare-gas cluster cation, Rg^, con- 
sists for the SO-free model of 3n valence-bond Slater determinants, \^k,p„^) 
(where k = 1, • • • , n and m = x,y, z) and represents states with the posi- 
tive charge localized in a valence pm-orbital of atom k. The corresponding 
3n X 3n hamiltonian matrix, 

is constructed as described in Ref. [Ij from diatomic potential energy curves 
for the electronic ground state of the neutral dimer, Rg2, and the electronic 
ground and three lowest excited states of the ionic dimer, Rg^. In our 
calculations, we have used semiempirical curves for neutral dimers [SI [7] and 
accurate ab initio curves for ionic diatoms [U |9] . 

This simple model is easily augmented [5] with SO coupling terms via a 
semi-empirical atoms-in-molecules scheme [10]. If the SO coupling is taken 
into account, the number of the basis set wave functions doubles, as there 
are two possible orientations of the spin of the electron removed from the 
valence shell of a particular atom, Sz = ±1/2, as well as the dimension of 
the electronic hamiltonian matrix. In addition, the matrix of Eq. [2] must be 
replaced by {6 denotes the Kronecker delta) [5J 

^k,pm,s,;k',p'^,si - -«fc,p™;fc',p^0s,;4 + ;p^,4 0^;^' ' W 

where 

f^f^a,lJk;k' = C{(l>k,p^,sjLkSk\(l)k',p',„,sO^ (4) 

^ is the SO coupling constant, and Lk and Sk are angular and spin operators 
for the k-th atom, respectively. The SO constants are independent inputs to 
the model and have been extracted here from experiments reporting on the 
SO splitting between the and ^^3/2 states of atomic monomers |llj . 

Hereafter, we denote this extended DIM model by DIM+SO. 

In the following text we use a simplified indices for the electronic wave 
function and hamiltonian matrix components, e.g., a = [k,pm,Sz] etc. 

'^In our work, the zero of energy is identified with the energy of fully dissociate state, 
Rg^ + (n — l)Rg , calculated with the SO interaction not included and the atomic con- 
tributions of Eq. [1] can be omitted. If, on the other hand, the SO coupling is considered, 
the atomic contribution corresponds either to the energy of Rg+('^P3/2) or Rg+(^Pi/2) 
measured from the SO-free Rg"*" level. 
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3 Non-radiative dynamics 



The semi-classical dynamical method (classical nuclei and quantum elec- 
trons) we use in our work for the non-radiative stage of our calculations, 
the MFQ-AMP/S method of Ref. [2j, combines a) the well known Ehrenfest 
mean- field approach [T2] , detailed for the rare-gas cluster cations in fl] , with 
b) the inclusion of quantum decoherence as introduced in Ref. \2l- 



3.1 Mean-field method 

The equations of motion for a system of classical nuclei surrounded by a 
cloud of electrons can be written within the mean-field approximation as 
coupled classical Hamilton equations for the nuclei 

Qi = — , (5) 

K = (V'|-^|V'> (6) 
and time dependent Schrodinger equation for the electrons 

.n'-f^m- (T) 

In Eqs. [5]-[7j and pi denote respectively generalized nuclear coordinates 
and momenta, H denotes the electronic hamiltonian, which depends para- 
metrically on the nuclear coordinates, and is a time dependent wave 
function representing the current electronic state. Small latin indices are 
used to label nuclear degrees of freedom and range between 1 through 3n. 

Within the DIM-I-SO approach, the electronic wave function, \ip), can be 
expanded using basis set wave functions of Sec. [2| \^a), also parametrically 
dependent on nuclear coordinates and, consequently, on time as well, 

a 

with a introduced above, a = [k,pm,Sz]- The electronic hamiltonian can 
also be expressed in an expanded form 

R = J2Hp^\^p){^,\, (9) 

/3,7 

where Hj^^ = Sp^Hfj^xSxj (with H^^x being the DIM-I-SO hamiltonian matrix 
given by Eq. |3j for simplicity we omit the (SO) upper index), and 5^^ = 
($Q,|$^) are overlap matrix elements. Note that matrix i?^^ is equal to the 
DIM-l-SO hamiltonian matrix, H^-y, if the overlaps are neglected (5q,^ = 
for a (3) and wavefunctions \^a) are normalized {Saa = 1)0 

^This is a usual and sufficiently accurate approximation adopted in all DIM models as 
yet developed for the rare-gas ionic clusters. 
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After inserting the expanded forms of the electronic hamiltonian and 
time-dependent electronic wave function into Eq. [6| one obtains 



Pi = - 



(10) 



where D 



= (<I)q,|^^) are non-diabatic coupling coefficients, and asterisks 



denote complex conjugation. The overlaps and non-diabatic couplings are 
usually neglected in DIM approaches and, consequently, Eq. 10 can be 
further simplified by setting Sap 
Kronecker delta), 



and D 



a,p "^"^ 



0, (with 5ap being the 

Pi = - } ^a-attp^^, (11) 

where, after neglecting the overlaps, Hap is replaced with Hap- 

Further simplification of Eq. 11 is possible if coefficients Cq, and matrix 
elements H^p, which are in general complex, are rewritten using their real 

and imaginary parts, Oq 



at^ + iJt^ and H^p 



Pi 



E 

a,p 



,(re)^(re) 



(im)„(im) 

/3 



dH 



(re) 



dqi 



+ [a. 



(re)^(im) _ „(im)^(re) 



'P 



'P 



dH, 



(im)- 

aP 



dqi 
(12)- 



The imaginary part of the electronic hamiltonian matrix is non-zero only if 
the SO coupling is included. If it is done using the Cohen-Schneider, atoms- 
in-molecules scheme p2| , all the imaginary terms are constant as they do not 



ap 



depend on the nuclear positions, and, consequently, = 0. The second 

term on the right-hand-side of Eq. [12] thus vanishes and the equation can 
be written in the final form 



Pi 



(re)„(je) ^^(im)„gm) 



"P 



'P 



a,P 



dH, 



(re) 



dqi 



(13) 



Similarly, the electronic Schrodinger equation, Eq. [7j can be rewritten 
after inserting the expansion of Eq. [8] to 

d\^p) 



E«/3H|^/3>> 

/3 



and after multiplying by {^a\ from the left to 

ih^Sapap + i^YDl^papqj = ^ H^sap. 



(14) 



(15) 



^ It directly follows from hermicity of the electronic hamiltonian matrix that its real 
part, H^l^\ is symmetric and the imaginary part, -ff^™', is antisymmetric. After using 
this property, we obtain immediately Eq. [12 
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A significant simplification of Eq. 15 is further possible if the overlap matrix 



considered diabatic, D^'^l = 0, 



is replaced with the Kronecker delta, Sa/s = Sap, and basis set \^a) is 

,(0 

a/3 

ifiaa = ^Hapap. (16) 



Eq. [16] must be treated with care, however, namely due to rapid oscil- 
lations occurring in the electronic wave function and, consequently, also in 
expansion coefficients aa- A special scheme has been developed for tack- 
ling this problem in the previous work. Since it is of technical rather than 
methodological importance, it is not discussed here and the reader is directed 
to Ref. [T| for details. 

3.2 Inclusion of quantum decoherence 

As shown elsewhere [2j, quantum decoherence is important, particularly for 
the heavy rare gases, krypton and xenon. It is introduced into the mean-field 
approach by periodically quenching the electronic wave function]^ We de- 
note this extended dynamical approach by MFQ (Mean Field with Quench- 
ings). The quenching algorithm comprises basically two steps. Firstly, the 
probabilities for collapsing the current electronic wavefunction into one of 
adiabatic states is calculated and a wave function collapse is proposed ac- 
cording to these probabilities. Secondly, in case the proposed collapse has 
been accepted, the kinetic energy of nuclei is adjusted so that the total 
energy of the system remains unchanged. The proposed jump can be, in 
general, rejected in both steps of the present algorithm and, in that case, 
the system resumes the coherent evolution until the next hop attempt]^ 

Several quenchings schemes have been developed previously [5]. In this 
work we use computationally cheap, but several times successfully tested 
MFQ-AMP/S algorithm. The procedure starts with calculating the adia- 
batic amplitudes of the current electronic wave function (hence the acronym 
AMP). More specifically, the normalized probability for collapsing the cur- 
rent electronic state, tp, to a particular adiabatic state, 0^, is calculated 
as 

(17) 

where represents the diagonal element of the electronic density matrix 
{Pfiu = c^cl and c^^ are amplitudes of the current electronic wave function. 



''Proper settings of the quenching period was thoroughly discussed in 2 . In this work 
we use quenching period fquonch = 100 fs. 

^For the algorithm used in this work, AMP (see below), the proposed electronic jump 
is always accepted in the first step and rejection can occur only during the second step, 
namely, if there is not enough kinetic energy to cover expenses of an upward electronic 
jump. 
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ip, expanded in the adiabatic basis set, tp = ^^c^ct)^). After the electronic 
jump is complete, the kinetic energy of nuclei is adjusted so that the total en- 
ergy of the system is conserved. In the MFQ-AMP/S model, this is achieved 
by scaling (hence the third acronym, S) nuclear velocities, as rationalized in 

4 Radiative dynamics 

After the non-radiative dynamics is stopped at time tDD) each trajectory 
is evaluated as an ensemble undergoing first-order decay due to radiative 
transitions in the electronic subsystem. In principle, many decay processes 
may occur in such an ensemble, both parallel and serial, which may lead to 
a complex system of coupled first-order equations governing the time evo- 
lution of this ensemble. In principle, such equations can be derived and 
solved. Nevertheless, since in our case a) transitions are expected only from 
the upper family of states of the charged fragment, an excited state resulting 
for the particular trajectory from the non-radiative dynamics at toD; to the 
lower family of states and b) the fragments undergo, after the radiative tran- 
sition, a rapid non-radiative decay, the radiative processes can be assumed 
parallel and a simplified set of decay equations can be used. In particular, 
if a population of n/o identical initial states from the upper family of states 
(e.g., all being state /) is considered for a particular trajectory and assumed 
to decay to the lower family of states (J), the corresponding population 
numbers will change with time At = t — too according to (dot denotes the 
time derivative) 

n/ = - ^ Turn, hj = Turn, (18) 
J 

with initial conditions n/(At = 0) = n/o (=1 for one particular trajectory) 
and nj{At = 0) = 0. It is easy to find, that the only solution to these 
equations is given by Eq. (2) of the letter, namely, 

nj(At) = n/o(l - e-r^*)r,j/r, ni{At) = nioe-^^\ (19) 

where F = X^x^i ^IK and At w t since t ^ tDD- 

Note also that nj{At)\niQ=i gives the probability that, at time At, the 
system will be found in state J, and n/(At)|„^Q=i is the probability of sur- 
viving the system in excited state /. The evaluation of fragments at time 
At consists then in a cycle repeated for all trajectories and comprising the 
following steps: 

1. identify the fragmentation channel corresponding to the particular tra- 
jectory, 

2. subtract from the total number of trajectories leading to the same 
fragmentation channel at tDD value of 1 — nj{At)\njQ=i, 
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3. identify the fragmentation channel for each state J (this can be done 
either by running additional non-radiative dynamical simulation or by 
simple energetic considerations), 

4. add nj{At)\njQ=i to the number of trajectories leading at ^dd to the 
same fragmentation channel. 

After this cycle is complete, one gets updated abundances of fragments as 
should be detected at time At f« t of radiative decay. 

The decay rates of Eq. [18] are calculated from a standard formula for 
spontaneous radiation (Eq. 1 of the letter), 

T/j = ^ 3 (El - Ejf Iml' , (20) 

ill eon 

where the transition dipole moment is obtained for a particular charged 
fragment geometry, R, within the point-charge approximation |13j . 

n z +1/2 

Mz.(R)«eEE E OgLl^^- (21) 

k=lm=x s^=-l/2 

The first sum of Eq, [21] runs over all atoms in the charged fragment and 

'^^kpmSz ^^'^ ^ipLsz ^'^^ amplitudes of adiabatic states / and J, respectively, 
expressed in the DIM+SO basis set introduced in Sec. [2] Alike in our 
earlier work, the point-charge approximation has been further improved in 
the present work by including damped polarization effects |14] consisting in 
a replacement 

Rfc Rfc E "efr(^jfc)-^T^, (22) 

where Rj^ = Rj — Rj, Rik = |Rifc|, and alg{Rik) is a damped effective 
polarizibility expressed in atomic units [14J . 
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